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We establish a general equivalence between van der Waals interaction energies within the formal- 
ism of the non-local van der Waals functional of the density functional theory and within the formal- 
ism of the field approach based on the secular determinants of the electromagnetic field modes. We 
then compare the two methods explicitly in the case of a planparallel geometry with a continuously 
varying dielectric response function and show that their respective numerical implementations are 
not equivalent. This allows us to discuss the merits of the two approaches and possible advantages 
of either method in a simple model calculation. 
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I. INTRODUCTION 

Van der Waals (vdW) interactions are ubiquitous and 
their effects are quintessential for understanding various 
properties of physical, chemical and biological systems 
In particular the planparallel layer geometry has 
played a prominent role in our understanding of phos- 
pholipid 0] and polymer assemblies 0, H| , as well as in 
inorganic systems such as the intergranular films in sili- 
con nitride structural ceramics [H or interfaces and grain 
boundaries in perovskite based electronic ceramics [6|]. 
Understanding molecular interactions in these systems is 
an important step in controlling the assembly process. 
Though interactions in these assemblies are due to many 
different specific properties, vdW interactions are a com- 
mon underlying and unifying feature. 

The vdW interaction energy (or free energy if temper- 
ature effects are of relevance) can be derived in many 
different ways [l|, 0] which are in principle equivalent but 
vary in terms of their implementability for specific prob- 
lems and geometries. Though all these approaches differ 
in various ways, they are all based on considerations of 
the electromagnetic field equations and the changes in 
the configurations of these fields when bodies are brought 
into close proximity. 

A different approach to vdW interactions stemms from 
the application of density functional theory (DFT) Q 
that has been very successful in the study of individ- 
ual molecules as well as dense solids. In this frame- 
work the vdW interactions belong to the non-local ef- 
fects that need to be considered when applying DFT to 
sparse matter [9ij. It is now well recognized that 
long-range effects in the DFT are not accounted for by 
the local or semilocal approximations of the exchange- 
correlation functional, though the vdW interaction is 
indeed a long-range correlation effect that is neverthe- 
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less unrelated to exchange correlations. A large body of 
work (see d, [13, [0, EH and references therein) has thus 
been devoted to a seamless extension of the DFT energy 
functionals that would properly capture the saturation 
of the vdW interaction at small separation and its con- 
tinuous approach to the standard DFT short-range in- 
teraction energy. In this context the vdW interactions in 
the planparallel slab geometry are derived via the matter 
response functions and successfully capture the seamless 
transition of long-range vdW interactions to those im- 
portant at near contact conditions [l2T |. 

In this work we will show that formulations of vdW 
interactions based on the field point of view [l3[ and on 
the matter point of view 0] give the same interaction en- 
ergy in the planparallel slab geometry with continuously 
varying dielectric response. This is just a corrolary on 
the general equivalence between the field and the matter 
approaches to the vdW interactions. Apart from that we 
will also show that the two approaches are however not 
equivalent when it comes to numerical implementation of 
the calculation of the vdW interactions, where the field 
approach seems to be numerically advantageous. 



II. DESCRIPTION OF THE PROBLEM 

We first present the expression for the vdW interaction 
energy derived within the DFT approach and the field 
approach and later show their equivalence. In what fol- 
lows we will limit ourselves to the nonretarded limit of a 
continuously (with a continuous derivative) varying local 
dielectric profile in the coordinate perpendicular to the 
two apposed planar surfaces at zero temperature. There- 
fore we can write for the dielectric response function 



e(r, r'; i£) = S(z-z') 



d 2 Q 

(2tt)2 



e(Q;z;i£,) e 



zQ(p-p') 



(1) 



where Q is the wavevector perpendicular to the axis z, 
p is the 2D radius vector in the plane perpendicular to 
z axis, z£ is the imaginary frequency with the dielectric 
function of imaginary frequencies being real and decreas- 
ing with increasing frequency [l[. In what follows we 
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with 5zi the slab thickness, i.e. the separ ation between 
the i—th and (i+l)—th layers. Here pi = \J Q 2 + e^/c 2 , 
and the quantities describing the relative dielectric miss- 
match at each layer i arc defined as 



for TM and 



Pi+l^i + Pi^i+1 



FIG. 1: A schematic diagram of a typical dielectric profile 
(arbitrary scale) and the interval [0, L] as required by the 
DFT approach. Near both ends of the interval the dielectric 
response should be constant for about at least oc 1/Q. 



shall use e(z) as a shorthand notation for e(Q;z;i£). 
Within the DFT formalism proposed by Rydberg et al. 
[ll[ the vdW interaction energy can be obtained as 



<PQ D M (i£;L) 



(2n) 2 



(2) 



where the expression Du(i^L) referrs to the two pla- 
nar interfaces at separation L, while D$(iQ referrs to 
a reference system (usually empty space). In the zero- 
temperature limit, the computation of Z?a/(i£;L) and 

Dftf(i£) can be reduced 0, [ll[ to the solution of the 
modified scalar Laplace equation 



-n . e -/ 
tp +—tp 

e 



Q 2 ¥> 



0. 



(3) 



if one identifies 
D M (i£;L) = (0'(O)) 



with boundary conditions 

The origin and the distance L need to be chosen such that 
the dielectric constant does not vary sufficiently close to 
both edges of the interval (within a distance oc 1/Q), as 
is shown in Figure [TJ 

In the field approach the calculation of the vdW in- 
teraction energy can be reduced to an algebra of 2 x 2 
matrices fllij . The interaction energy still retains the 
form 



m-i | 



d 2 Q ^D F (it,;L) 



(2^) 2 Df(*0 



(5) 



but here the interpretation of D(i£; L) is different. In 
fact one derives that Dp(i^; L) = M u is the (1, 1) matrix 
element of the product of matrices 



M 



■ Ti + \Di + iTiDi 



(6) 



where the index i runs over the whole interval of properly 
discretized z axis, defined as 



Ti = 



1 

exp(— 2piSzi 



A: = 



1 

-A,; 



-A 4 
1 



(7) 



for TE modes, where we assume that the magnetic per- 
meability equals 1 in the whole spatial domain. Tj matrix 
corresponds to the phase shift of the left and right trav- 
elling waves when moving across the i-th slab, and 
matrix links the fields at the boundaries between slabs. 
The matrices T and D are strictly real. 

Off hand it is not clear whether the two expressions, 
i.e. (H]) and (J3J), for the vdW interaction energy are the 
same, since the definitions of D(i£;L) in both cases are 
not obviously related. In what follows we will first of all 
show that the two expressions are not only related but 
are, modulo some unimportant scaling factors, in fact 
exactly the same. 



III. PROOF OF EQUIVALENCE 

We now investigate whether the two approaches to the 
vdW interaction energy give formally the same result. 
We start with equation ([3]) from the DFT approach and 
note that the solution of any equation of the form 



y" + (3(z)y'-Q 2 (z)y = 



#0) = o, m) = i. 

can be written as 



y = yx + 2/2, 



(8) 



(9) 



where y% and yi satisfy a system of two coupled linear 
differential equations 
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Q(z)-j(z) j(z) 
j(z) -Q(z)-f(z) 
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with 



7(*) 



Q'(z) 
Q(z) 



(10) 



(11) 



In the DFT approach, as Q does not depend on x, we 
have 



70) = i(ln( e ))'. 



(12) 



We will now try and reformulate the matrix product 
© in the field approach so that it will lead to the same 
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equation (fTO)) as the DFT approach. The matrix ap- 
proach can be somewhat simplified by noting that the 
initial condition for the matrix M can be written as 



This is still not equal to eq. (jlO)) . The next transforma- 
tion necessary is 



M 



' 1 


" 




' 1 " 














® [ 1 



(13) 



Ipi =1pl, 1p2 

which now yields finally 



This is due to the fact that, at nonzero p, the product 
of matrices Tj for the homogeneous space extending to 
— oo will converge to the above matrix, and thus it can 
be taken as the initial condition for the matrix M at any 
coordinate zq, where for z < zq the dielectric constant 
is nonvarying. Such splitting of the matrix reduces the 
problem of the matrix product to a problem of matrices 
acting on a vector, in the sense that after i layers 
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Q~l{z) 

7(z) -Q-70) 



4>2 



(20) 



(21) 



Obviously now the two sets of equations (fTU|) and (JT7J) 
coincide, which is the main step towards demonstrating 
the equivalence. Note also that the transformation 

Xl (z) = exp( f (Q + j(z'))dz') iPi(z) = 



T i D i T i - 1 D i - 1 ...TiD 1 



(14) 




exp(Qz) iPi(z) (22) 



At the final layer we may identify ipi = Mn, ip2 = M21 
and therefore 



yields the system of equations 

Q 







"xT 




Df = ipi ■ 


(15) 


x' 2 _ 





- 7 (z) 



-j(z) 

Q + i{z) 





Xi 




. * 2 . 



(23) 



Let us now reduce the matrix product in the field for- 
mulation into a set of coupled differential equations for a 
continuously varying dielectric profile. Here we will treat 
only the nonretarded case, c — ► 00, where p = Q. If the 
layer thickess is small, both the matrix Ti as well as Di 
differ only slightly from the identity, such that 



According to equation (fT0|) , the sum H± — \i + Xi then 
satisfies the differential equation 



H'l - -H' ± - Q 2 H ± = 



(24) 





-2Q 



Szi, Di w 7+ 





-It 



-It 




Sz h (16) 



where 7, = Ai/Szi. This is a difference scheme that in 
the limit Szi — > leads to a continuous formulation of 
the equation (fT4")l . 



^2 





- 7 (z) " 




ipi 




- 7 (z) -2Q _ 




1p2 



and describes the TM magnetic field modes. The only 
difference w.r.t. equation ([3]) is the sign in front of the 
term containing the first derivative, but we will show that 
this sign difference does not affect the interaction energy 
computation. 

We thus need to show that the expressions dU) and 
(fT5)) for D F and D M as used in the energy integral in 
equations (J2|) and (J5J) give the same result. Let us rewrite 
the energy expression in the field approach as 



(17) 



In 



with the initial condition ipi(-oo) — 1, ^2(^00) = 0. 
Let us now transform the system of differential equa- 
tions fTT)) obtained via the matrix product calculation so 
that it will coincide with the system (fit)]) as derived from 
the DFT approach. Obviously eqs. JTUI) and (JTTJ) differ 
in the diagonal elements of the matrix. However, if we 
introduce the transformation 



^(z) = exp ( J (Q - j(z ))dz )ipi(z) 



Df{L) 
Df[L) 



= In 



ML) 



= In 



ML) 4 0) (Q) 
Mo) ^i\l) 



(25) 



exp( Q z ) ^ (z ) (is) 



we then derive the following system of equations 
$'2 





Q-~j{z) -7(2:) 




'M 




-j(z) -Q-j{z)_ 




1p2 



(19) 



Since the additional terms V'i(O), ip^\o) are equal to 
1 with the usual choice of initial conditions, due to the 
linearity of the problem, the above expression allows for 
a scaling of the solution without changing the result. 

It now helps to think about the DFT approach in terms 
of the split system of equations (fTU|) , where (p = y\ + y%. 
As the second component y2 is exponentially decreasing, 
when L is sufficiently large, we may write <p(L) w y\ 
and therefore yi(L) = 1. For the boundary condition 
(p(0) — to be fulfilled, we must require yi(0) = — 2/2(0). 
This is, of course, a different boundary condition than in 
the matrix approach, where V'i(O) — L ^2(0) = 0. Since, 
however, in a homogeneous region the second component 
decays exponentially fast, and if the homogeneous region 
extends sufficiently far into the z > region, the two 
solutions are largely equivalent when rescaled properly. 
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According to equation (|10|). the derivative that enters 
the free energy computation is then <p'(0) = 2Q ?/i(0). It 
then follows that 



In 



Dm(L) 



= In 



^(0) 



2Q y?\o) = 
2Q 2/i (0) 

Vi(L) yj 0) (0) 



= In 



In 



(26) 



which agrees with the expression f|25[) save for the trans- 
formations (fT5|) and (|20[) which take us from the variables 
ipi to m. 

To show that these transformations do not affect the 
result, let us now assume that the dielectric functions e(z) 
for the system at study and the reference system have 
the same limit at z — > oo, while at the same time they 
also both limit to a possibly different value at z — » -co. 
If the transformations (fTSl [20)) are then applied to the 
expression (|25j) . the effects of the transformations exactly 
cancel for the original and reference systems and we are 
left with the expression (|2"6"|) , concluding the complete 
demonstration of equivalence. 



IV. CONSEQUENCES AND DISCUSSION 

Though we have just shown that in principle the DFT 
and the field approaches are equivalent, there are subtle 
differences in the numerical implementation that allow 
for a sensible comparison of the two methods which we 
discuss below. 

First we would like to comment on the Lifshitz limit, 
where the dielectric profile consists of two semi-infinite 
halfspaces with one dielectric constant, separated by a 
finite slab of thickness d with a different dielectric con- 
stant. While the matrix product approach can deal with 
problems of this type directly, the continuous picture of 
either the field or the DFT approach encounters difficul- 
ties. The first issue is the definition of the derivative 



7 (z) = i(ln( £ ))' = ^ 



(27) 



at the dielectric boundary. While for an interface at z = 
and having a jump in the dielectric constant Ae we may 
write e' = Ae 5(z), it is the value of e in the denomi- 
nator of the above equation that is not well defined at 
the interface. Furthermore, equations (f2"4"|) and (J3J) put 
7 and therefore the S(z) function in front of the first or- 
der derivative term, making it impossible to establish the 
proper connection formulae without further physical in- 
sight. This insight is, however, already incorporated in 
the discrete version of the matrix product formulation in 
which the Lifshitz case is handled gracefully and emerges 
naturally. 

In order to demonstrate the concepts outlined in the 
previous section we introduce a test case of the dielectric 




FIG. 2: Plot of (In e) for the model dielectric response func- 
tion, equation (|28[1 . for the values of parameters parameters 
a = 1, a = 1/4 and A = 2. 




FIG. 3: Comparison of the solution of the two differential 
equations (|24|) and ((3| that only differ in the sign of the first 
order derivative term. In both cases, the result was divided 
by oc exp(Qz). The full line shows the positive sign as in 
equation ©, the dashed line the negative one as in equation 
(|24p . The dotted line shows the value of 1. 



profile at a certain imaginary frequency in the integral 
©, given by 



(hie)' = A 



exp 



(z - af 



exp 



{z + af 



2a 2 J 1 V 2ff2 

'(28) 

as shown in the figure [2] with parameters a — 1 , a = 1/4 
and A = 2, which are also the chosen parameters in all 
the calculations shown further. In the figure [3] we show 
the numerical solutions of the differential equations ([24]) 
and that differ in the sign of the first order deriva- 
tive term. All calculations are done for Q = 1. What 
is shown for each of the differential equations is the ra- 
tio of the solution starting with rather arbitrary initial 
conditions at large negative z (ideally z — > — oo) to the 
solution of the same differential equation with the same 
initial condition but for empty space, which itself behaves 
as oc exp(Qz). While the rescaled solutions of the two 
differential equations can be seen to deviate markedly 
from one another in the region of the varying dielectric 
constant as well as in the region of the order of oc 1/Q 
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FIG. 4: The solution obtained by solving the second order 
differential equation (full line, equation ((3} rescaled as in the 
figure [3]) directly, as opposed to the solution of the system of 
the coupled first order differential equations in equation (|I0[) . 
The dashed line shows the first component of the coupled 
system, the dashed-dotted line shows the second component 
plus one, with the circles representing a direct sum of the two 
solutions. The thin dotted line represents the value of 1. The 
two complete solutions are indistinguishable from each other. 

beyond it, when z — > oo they relax to the same value, 
and for the calculation of the free energy Q this value 
is the only important quantity. In figure [4] the compari- 
son of the direct solution of the differential equation ([3]) 
with the solution obtained by splitting it into the system 
of differential equations (fTTj|) is shown. Again, the re- 
sults are divided by the empty space solution cx exp(Qz). 
The sum of the two components of the coupled system 
of equations is indistinguishable from the direct solution 
to the original equation. What is apparent is that the 
solution for the first component of the coupled system 
of equations varies only in the region of the varying di- 
electric constant (see the figure [2]); beyond that region 
only the second component of the system is still relaxing 
to with the first component remaining constant. For 
the purposes of the calculation of the free energy, the use 
of a coupled system is therefore beneficial as the calcu- 



lation can be terminated immediately after entering the 
homogeneous region, seeing that the first component of 
the system carries all the information about the required 
asymptotic behaviour, whereas in the direct solution to 
the differential equation it is necessary to go deeper into 
the homogeneous region for the solution to relax. 

This shows that using the matrix product formulation 
based on the field approach has several advantages over 
the DFT approach in the numerical implementation of 
the calculation. The main benefit is that the domain of 
calculation necessary to determine the solution is only 
that part of the z-axis beyond which the dielectric re- 
sponse is constant, whereas in the continuous DFT ap- 
proach one needs to extend the system by the length 
cx l/Q in order to allow for the relaxation of the solution. 
Furthermore, the matrix product approach allows for di- 
rect treatment of discontinuous dielectric profiles which, 
while in principle accessible via the DFT approach, nev- 
ertheless require a somewhat special treatment. The ad- 
vantage of the field approach is particularly relevant for 
cases where there are several discontinuities or mixed 
continuous regions interspersed with discontinuities in 
the dielectric response function that can all be handled 
quite straightforwardly in the field approach. Finally, 
while we showed the equivalence of the two approaches 
only for the nonretared case of the Van der Waals inter- 
actions, valid at small separations between interacting 
systems, the field approach can be generalized directly 
and straightforwardly to calculate the retarded van der 
Waals - dispersion interactions as well. As far as we know 
the DFT approach has not yet been generalized into the 
retarded van der Waals region. 
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